function Mesh = Mesh_struct_2d(PDE, Nx, Ny, basis_type)
Mesh.left = PDE.domain(1);
Mesh.right = PDE.domain(2);
Mesh.bottom = PDE.domain(3);
Mesh.top = PDE.domain(4);
Mesh.Nx = Nx;
Mesh.Ny = Ny;
Mesh.hx = (Mesh.right-Mesh.left)/Nx;
Mesh.hy = (Mesh.top-Mesh.bottom)/Ny;
Mesh.basis_type = basis_type;
[Mesh.node, Mesh.elem] = process_node_elem_2d(Mesh.left, Mesh.right, Mesh.bottom, Mesh.top, Nx, Ny, basis_type);
Mesh.Np = (Nx+1)*(Ny+1);
Mesh.Ne = size(Mesh.elem,1);
Mesh.Nb = size(Mesh.node,1);
Mesh.Nlb = size(Mesh.elem,2);
Mesh.absJ = Mesh.hx*Mesh.hy;
end

%% process_node_elem_2d
function [node, elem] = process_node_elem_2d(left, right, bottom, top, Nx, Ny, basis_type)
[X, Y] = meshgrid(linspace(left, right, Nx+1), linspace(bottom, top, Ny+1));
node = [X(:), Y(:)];
Nb = size(node,1);
k = (1:Nb-Ny-1)';
cut = (Ny+1)*(1:Nx);
k(cut) = [];
elem = [k, k+Ny+1, k+1;
    k+1, k+Ny+1, k+Ny+2];
elem = sortrows(elem);
switch basis_type
    case "P1b"
        new_node = zeros(size(elem,1),2);
        for k = 1:size(elem,1)
            V = [node(elem(k,:),:)];
            new_node(k,:) = polygon_centroid(V);
        end
        node = [node; new_node];
        new_elem_col = (Nb+1:Nb+size(elem,1))';
        elem = [elem, new_elem_col];
    case "P2"
        [X, Y] = meshgrid(linspace(left, right, 2*Nx+1), linspace(bottom, top, 2*Ny+1));
        new_node = [X(:), Y(:)];
        cut = zeros(Ny+1,Nx+1);
        k = 1;
        for i = 0:2:2*Nx
            cut(:,k) = i*(2*Ny+1)+1:2:(i+1)*(2*Ny+1);
            k = k+1;
        end
        new_node(cut(:),:) = [];
        node = [node; new_node];
        t = Nb+1;
        c = 1;
        for k = 1:2:size(elem,1)
            i = mod(c-1,Ny);
            new_elem_col(k:k+1,:) = [t+Ny+i, t+Ny+1+i, t;
                t+Ny+1+i, t+3*Ny+1, t+Ny+2+i];
            if mod(c,Ny) == 0
                t = t+2*Ny+2;
            else
                t = t+1;
            end
            c = c+1;
        end
        elem = [elem, new_elem_col];
end
end

%% polygon_centroid
function center = polygon_centroid(V)
n = size(V,1);
A = 0;
center_x = 0;
center_y = 0;
for i = 1:n
    j = mod(i,n) + 1;
    xi = V(i,1);
    yi = V(i,2);
    xj = V(j,1);
    yj = V(j,2);
    factor = xi * yj - xj * yi;
    A = A + factor;
    center_x = center_x + (xi + xj) * factor;
    center_y = center_y + (yi + yj) * factor;
end
A = 0.5 * A;
center_x = (1 / (6 * A)) * center_x;
center_y = (1 / (6 * A)) * center_y;
center = [center_x, center_y];
end